#include "uhm.hxx"

#define N 1

using namespace uhm;

typedef ID_Node_H_<N>          Id_Nd;
typedef ID_Element_<N>         Id_El;

typedef LU_nopiv_Linal_        Uhm;                                                                         
//typedef LU_piv_Linal_           Uhm;
//typedef LU_incpiv_Linal_        Uhm;   

typedef Direct_Solver_<Id_Nd,Id_El,Uhm> Solver;

int main (int argc, char **argv) {

  initialize(argc, argv);

  // ** Environments
  Int_
    type = Mat_Base_::DOUBLE,
    uplo = Mat_Base_::FULL,
    sym  = Mat_Base_::UNSYM,
    rhs  = 1,
    blk  = 256;

  Int_
    kind = Mesh_Base_::NODE_DEFAULT;

  // ** Solver initialization
  Solver solver(type, uplo, sym, rhs, blk);

  // ** Set up caching level
  option_begin();
  Int_ caching;
  get_option_int("-caching", "Setting caching level", __FILE__,
                 0, caching);
  solver.set_caching(caching);
  Int_ n_iter;
  get_option_int("-iter", "Number of solving iteration", __FILE__,
                 2, n_iter);
    std::string fin;
  get_option_string("-in", "Reading a file", __FILE__,
                    "../../../uhmfile/toy_1.uhm", fin);
  Int_ time;
  get_option_int("-time", "Time tracking level 1) User 2) Front End 3) Internal",
                 __FILE__,UHM_TIME_LEVEL_FRONT,time);
  option_end();

  time_begin(time);

  for (Int_ i=0;i<n_iter;++i) {
    std::cout << "Caching level = " << caching << std::endl;
    solver.set_caching(caching++);
    
    std::cout << "Opening file ... " << fin << std::endl;
    std::ifstream in;
    in.open(fin.c_str());
    stream2mesh<Solver::_Mesh_>(in, solver.mesh());
    in.close();
    
    UHM_ERROR((mesh2cache<Solver::_Mesh_,Solver::_Cache_>(solver.mesh(), solver.work())),
              ">> Error in mesh2cache");

    if (i==0) {
      solver.work().set_matrices(type, rhs, blk);
      solver.work().create_matrices();
      solver.work().random_matrices();
    }

    UHM_ERROR(solver.create(),    ">> Error in Creating");
    
    time_in("updirect:: Decomposing");
    UHM_ERROR(solver.updecompose(), ">> Error in Decomposing");
    Double_ t = time_out();

    Double_ flop, recorded;
    solver.get_flop(flop, recorded);
    std::cout << std::endl << std::endl;

    std::cout << "FLOP = " << (flop/1.0e9)
              << ", Updated  = " << (recorded/1.0e9) << "  GFLOP"
              << ", Updating Percentage = " << (recorded/flop) << "%\n\n";

    std::cout << "FLOPS= " << (flop/1.0e9)/t
              << ", Modified = " << ((flop-recorded)/1.0e9)/t << "  GFLOPS\n\n";

    Double_ storage, max_storage;
    storage_counter(storage, max_storage);
    std::cout << "Storage = " << (storage/1.0e6) 
              << ", " << (max_storage/1.0e6) << "  MBytes\n";

    UHM_ERROR(solver.clear(),     ">> Error in Clearing");

    std::cout << i << " Iteration Done\n";
  }

  // ** Flush cache
  UHM_ERROR(solver.flush(),     ">> Error in Flushing");
  
  time_end();

  finalize();

  return true;
}


